package gis

/**
b, l, h : b: lat , l:lng


long double GetX(double RadLat, double eSquare, double a) {
        long double X;
        long double A0;
        long double A2;
        long double A4;
        long double A6;
        long double A8; ///*  from hsq
        A0 = 1 + 3.0 / 4.0 * eSquare + 45.0 / 64.0 * eSquare * eSquare
            + 350.0 / 512.0 * eSquare * eSquare * eSquare
            + 11025.0 / 16384.0 * eSquare * eSquare * eSquare * eSquare;
        A2 =
            -1.0 / 2.0
            * (3.0 / 4.0 * eSquare + 60.0 / 64.0 * eSquare * eSquare
                + 525.0 / 512.0 * eSquare * eSquare * eSquare
                + 17640.0 / 16384.0 * eSquare * eSquare * eSquare
                * eSquare);
        A4 = 1.0 / 4.0
            * (15.0 / 64.0 * eSquare * eSquare
                + 210.0 / 512.0 * eSquare * eSquare * eSquare
                + 8820.0 / 16384.0 * eSquare * eSquare * eSquare * eSquare);
        A6 = -1.0 / 6.0
            * (35.0 / 512.0 * eSquare * eSquare * eSquare
                + 2520.0 / 16384.0 * eSquare * eSquare * eSquare * eSquare);
        A8 = 1.0 / 8.0 * 315.0 / 16384.0 * eSquare * eSquare * eSquare * eSquare;
        X = a * (1 - eSquare)
            * (A0 * RadLat + A2 * sin(2 * RadLat) + A4 * sin(4 * RadLat)
                + A6 * sin(6 * RadLat) + A8 * sin(8 * RadLat));
        return X;
    }


    bool bl2Gaussxy(double RadLat, double RadLon, double* x, double* y,
        double RadLon0, double eSquare, double a, double heightchange) {
        long double m0, t, etaSquare, RadDeltaLon;
        long double N, X0;
        long double ePrimeSquare;
        long double aa, da, db; //ykoky2k
        RadDeltaLon = RadLon - RadLon0;
        ePrimeSquare = eSquare / (1 - eSquare);
        N = (a / sqrt(1 - eSquare * sin(RadLat) * sin(RadLat))) + heightchange;
        aa = N * sqrt(1 - eSquare * sin(RadLat) * sin(RadLat));
        da = aa - a;
        db = eSquare * sin(RadLat) * cos(RadLat)
            * (1 - eSquare * sin(RadLat) * cos(RadLat)) * da
            / (a * (1 - eSquare));
        RadLat = RadLat + db;
        a = aa;
        m0 = RadDeltaLon * cos(RadLat);
        t = tan(RadLat);
        etaSquare = ePrimeSquare * cos(RadLat) * cos(RadLat);
        X0 = GetX(RadLat, eSquare, aa);
        *x = X0 + 1 / 2.0 * N * t * pow(m0, 2)
            + 1 / 24.0
            * (5.0 - pow(t, 2) + 9 * etaSquare + 4 * pow(etaSquare, 2))
            * N * t * pow(m0, 4)
            + 1 / 720.0
            * (61.0 - 58 * pow(t, 2) + pow(t, 4) + 270 * etaSquare
                - 330 * etaSquare * pow(t, 2)) * N * t * pow(m0, 6)
            + 1 / 40320 * N * t * pow(m0, 8)
            * (1385 - 3111 * pow(t, 2) + 543 * pow(t, 4) - pow(t, 6));
        *y = 500000.0 + N * m0
            + 1 / 6.0 * (1.0 - pow(t, 2) + etaSquare) * N * pow(m0, 3)
            + 1 / 120.0
            * (5.0 - 18.0 * pow(t, 2) + pow(t, 4) + 14.0 * etaSquare
                - 58.0 * pow(t, 2) * etaSquare) * N * pow(m0, 5)
            + 1 / 5040.0 * pow(m0, 7)
            * (61 - 479 * pow(t, 2) + 179 * pow(t, 4) - pow(t, 6));
        return true;
    }

*/

import (
	"fmt"
	"math"
	"strconv"
	"strings"
)

const (
	WGS84_R  = 6.3781370e6           //长半轴(赤道半径(单位m))
	WGS84_FE = 1.0 / 298.257223563   // 地球扁率   /* earth flattening (WGS84) */
	WGS84_E2 = 0.00669437999013      // 第一偏心率平方
	WGS84_L  = 2 * math.Pi * WGS84_R // 赤道周长
)

const (
	_LNG_1KM    = 0.0095596713       // 经度1KM 的度数 (lng:75, lat:20)  CalcuPerM(lat) * 1000
	_LNG_1KM_55 = 0.0156617          // 经度1KM 的度数 (lng:135, lat:55)
	_LNG_1KM_30 = 0.0103728517984788 // 30纬度 1KM 的度数
	_LNG_1KM_31 = 0.010480045899962  // 31纬度 1KM 的度数
	//  _LAT_1KM = 0.0089831531890287     // 纬度1KM 的度数
	_LAT_1KM = 0.009009009009009 //  1度= 111KM 1/111
)

// 角度 -> 弧度
func D2R(d float64) float64 {
	return d * math.Pi / 180
}

// 弧度 -> 角度
func R2D(r float64) float64 {
	return r * 180 / math.Pi
}

// 计算纬度上
// 纬度为角度
//
//	返回该纬度上 每一度的周长 M 360 / 周长
func CalcuPerM(latD float64) float64 {
	latR := D2R(latD)
	r := WGS84_R * math.Cos(latR) // 该圆的半径
	l := 2 * math.Pi * r          // 周长
	return l / 360
}

func CalcuPerMV1(latD float64) float64 {
	latR := D2R(latD)
	l := WGS84_L * math.Cos(latR)
	return l / 360
}

// 3007.14645905  -> 30.11910765
func DDMM2Degree(val float64) float64 {
	r := val / 100
	v := math.Trunc(r) // 30
	r = v + (r-v)*10/6 // 30 + 0.0714 * 10 / 6.0
	return r
}

// // 30.11910765 -> 3007.14645905
func Deg2DDMM(deg float64) string {
	//d := math.Trunc(deg) * 100
	//m := deg - math.Trunc(deg)
	//m = m * 60
	//return fmt.Sprintf("%.8f", d+m)

	d := math.Trunc(deg)
	m := deg - math.Trunc(deg)
	m = m * 60
	if m < 10 {
		return fmt.Sprintf("%.2d0%.8f", int(d), m)
	} else {
		return fmt.Sprintf("%.2d%.8f", int(d), m)
	}
}

/*
*

	经度dddmm.mmmm（度分）格式(前面的0也将被传输)
	不支持负数, 负数先abs
*/
func Deg2DDDMM(deg float64) string {
	d := math.Trunc(deg)
	m := deg - math.Trunc(deg)
	m = m * 60
	if m < 10 {
		return fmt.Sprintf("%.3d0%.8f", int(d), m)
	} else {
		return fmt.Sprintf("%.3d%.8f", int(d), m)
	}

}

func Deg2DDMMSS_Lat(deg float64) string {
	return Deg2DDMMSS(deg, "S", "N")
}

func Deg2DDMMSS_Lng(deg float64) string {
	return Deg2DDMMSS(deg, "W", "E")
}

// 30.4334707->030°26′00.49452″N
func Deg2DDMMSS(deg float64, negative, positive string) string {
	suffix := ""
	if deg < 0 {
		suffix = negative
		deg = math.Abs(deg)
	} else {
		suffix = positive
	}
	d := math.Trunc(deg)
	m := deg - math.Trunc(deg)
	m = m * 60
	m0 := math.Trunc(m)

	s := m - m0
	s = s * 60

	s_str := fmt.Sprintf("%.5f", s)
	if len(s_str) < 8 {
		s_str = "0" + s_str
	}
	return fmt.Sprintf("%.3d°%.2d′%s″%s", int(d), int(m0), s_str, suffix)
}

// 030°26′00.49452″N -> 30.4334707
//
//	114°13′28.37346″E    114.22454818333334
//
// 30:26:00.51945  -> 30.433477625000002
//
//	114:13:28.36367       114.22454546388889
func DDMMSS2Deg(deg string) float64 {
	var d float64
	var m float64
	var s float64

	var dd string
	var mm string
	var ss string
	//int64, err := strconv.ParseInt(string, 10, 64)
	if strings.Contains(deg, "N") || strings.Contains(deg, "E") {
		dd = strings.Split(deg, "°")[0]
		d, _ = strconv.ParseFloat(dd, 64)
		mm = strings.Split(strings.Split(deg, "°")[1], "′")[0]
		m, _ = strconv.ParseFloat(mm, 64)
		ss = strings.Split(strings.Split(strings.Split(deg, "°")[1], "′")[1], "″")[0]
		s, _ = strconv.ParseFloat(ss, 64)
	}

	if strings.Contains(deg, ":") {
		dd = strings.Split(deg, ":")[0]
		d, _ = strconv.ParseFloat(dd, 64)
		mm = strings.Split(deg, ":")[1]
		m, _ = strconv.ParseFloat(mm, 64)
		ss = strings.Split(deg, ":")[2]
		s, _ = strconv.ParseFloat(ss, 64)
	}
	//fmt.Println("度",d)
	//fmt.Println("分",m)
	//fmt.Println("秒",s)
	return d + m/60 + s/60/60
}

func CRC(buf []byte) byte {
	if len(buf) == 0 {
		return 0
	}

	r := buf[0]

	if len(buf) > 1 {
		for _, v := range buf[1:] {
			r ^= v
		}
	}
	return r
}

// / 计算格网中间点
func CalcuGridCenter(lng, lat, gridlkm float64) (lngr float64, latr float64) {
	grdlng := _LNG_1KM * gridlkm
	grdlat := _LAT_1KM * gridlkm

	lngN := int(math.Trunc(lng / grdlng))
	latN := int(math.Trunc(lat / grdlat))

	lngr = float64(lngN)*grdlng + grdlng/2.0
	latr = float64(latN)*grdlat + grdlat/2.0
	return
}

// / 计算格网格编号
func CalcuGridNo(lng, lat, gridlkm float64) int64 {
	grdlng := _LNG_1KM * gridlkm
	grdlat := _LAT_1KM * gridlkm

	lngN := int(math.Trunc(lng / grdlng))
	latN := int(math.Trunc(lat / grdlat))
	return int64(lngN<<32 + latN)
}

// 计算网格中心点
func CalcuGridCenterByNo(gridno int64, gridlkm float64) (lngr float64, latr float64) {

	grdlng := _LNG_1KM * gridlkm
	grdlat := _LAT_1KM * gridlkm

	lngN := int(gridno >> 32)
	latN := int(gridno << 32 >> 32) // 去掉高位的32

	lngr = float64(lngN)*grdlng + grdlng/2.0
	latr = float64(latN)*grdlat + grdlat/2.0
	return
}

// / <summary>
// /  计算两点之间的距离(和GDMap计算北京距离一样 (南极，非洲，北京),
// /    http://lbs.amap.com/api/javascript-api/example/calcutation/calculate-distance-between-two-markers/
// /  return M
// /    http://www.cnblogs.com/zhoug2020/p/3950933.html
// / </summary>
func CalcuDistance(lng1, lat1, lng2, lat2 float64) float64 {
	lat1 = D2R(lat1)
	lng1 = D2R(lng1)
	lat2 = D2R(lat2)
	lng2 = D2R(lng2)

	if lat1 < 0 {
		lat1 = math.Pi/2.0 + math.Abs(lat1) // south
	} else if lat1 > 0 {
		lat1 = math.Pi/2.0 - math.Abs(lat1) // north
	}
	if lng1 < 0 {
		lng1 = math.Pi/2.0 - math.Abs(lng1) // west
	}

	if lat2 < 0 {
		lat2 = math.Pi/2.0 + math.Abs(lat2) // south
	} else if lat2 > 0 {
		lat2 = math.Pi/2.0 - math.Abs(lat2) // north
	}

	if lng2 < 0 {
		lng2 = math.Pi/2.0 - math.Abs(lng2) // west
	}

	x1 := WGS84_R * math.Cos(lng1) * math.Sin(lat1)
	y1 := WGS84_R * math.Sin(lng1) * math.Sin(lat1)
	z1 := WGS84_R * math.Cos(lat1)

	x2 := WGS84_R * math.Cos(lng2) * math.Sin(lat2)
	y2 := WGS84_R * math.Sin(lng2) * math.Sin(lat2)
	z2 := WGS84_R * math.Cos(lat2)

	d := math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2))

	//余弦定理求夹角
	theta := math.Acos((WGS84_R*WGS84_R + WGS84_R*WGS84_R - d*d) / (2 * WGS84_R * WGS84_R))

	return theta * WGS84_R
}

// / <summary>
// /  计算两点之间的距离(和GDMap计算北京距离一样 (南极，非洲，北京)
// /    http://lbs.amap.com/api/javascript-api/example/calcutation/calculate-distance-between-two-markers/
// /  return M
// /     http://www.jb51.net/article/61274.htm
// / </summary>
func CalcuDis(lng1, lat1, lng2, lat2 float64) float64 {
	if math.Abs(lng1-lng2) < 0.00000001 && math.Abs(lat1-lat2) < 0.00000001 {
		return 0
	}

	lat1 = D2R(lat1)
	lng1 = D2R(lng1)
	lat2 = D2R(lat2)
	lng2 = D2R(lng2)

	theta := lng2 - lng1
	dist := math.Acos(math.Sin(lat1)*math.Sin(lat2) + math.Cos(lat1)*math.Cos(lat2)*math.Cos(theta))

	return dist * WGS84_R
}

func CalcuDisXyz(x0, y0, z0, x1, y1, z1 float64) float64 {
	b0, l0, _ := EcefXyz2PosiLatLngH(x0, y0, z0)
	b1, l1, _ := EcefXyz2PosiLatLngH(x1, y1, z1)
	return CalcuDis(l0, b0, l1, b1)

}

func Round64(f float64, n int) float64 {
	p := math.Pow10(n)
	return math.Round(f*p) / p
}

func OffsetLatMetre(lat float64, n float64) float64 {
	v := n / 1000 * _LAT_1KM
	return lat + v
}

func OffsetLngMetre(lng float64, n float64) float64 {
	v := n / 1000 * _LNG_1KM_30
	return lng + v
}

func EcefXyz2PosiLatLngH(x, y, z float64) (lat, lng, h float64) {
	if x == 0 && y == 0 {
		return 0, 0, 0
	}
	lat, lng, h = Ecef2Pos(x, y, z)
	lat, lng = R2D(lat), R2D(lng)

	h = Round64(h, 5)
	return
}

func PosiLatLngH2EcefXyz(lat, lng, h float64) (x, y, z float64) {
	if lat == 0 && lng == 0 {
		return 0, 0, 0
	}
	return Pos2Ecef(lat, lng, h)
}

// /将 大地坐标系( lat, lon, h) (度) (WGS84) ->ECEF坐标系(X、Y、Z)转换成
func Pos2Ecef(lat, lon, h float64) (x, y, z float64) {

	lat_r := D2R(lat)
	lon_r := D2R(lon)
	N := WGS84_R / math.Sqrt(1-WGS84_E2*math.Sin(lat_r)*math.Sin(lat_r))
	x = (N + h) * math.Cos(lat_r) * math.Cos(lon_r)
	y = (N + h) * math.Cos(lat_r) * math.Sin(lon_r)
	z = (N*(1-WGS84_E2) + h) * math.Sin(lat_r)

	//var
	//	N:Double;
	//begin
	//b := D2R(b);
	//l := D2R(l);
	//
	//N := WGS84_a/sqrt(1-WGS84_e2*sin(b)*sin(b));
	//x := (N+h)*cos(b)*cos(l);
	//y := (N+h)*cos(b)*sin(l);
	//z := (N*(1-WGS84_e2)+h)*sin(b);

	return

}

///将 ECEF坐标系(X、Y、Z)转换成大地坐标系( lat, lon, h) (WGS84)
/*
X,y,z to b(r),l(r),h
  b, l, 为弧度
*/
func Ecef2Pos(x, y, z float64) (b, l, h float64) {

	e2 := WGS84_FE * (2.0 - WGS84_FE)
	r2 := (x * x) + (y * y) // dot(r,r,2)
	v := WGS84_R

	var sinp float64

	// 计算高 z
	_z := z
	zk := 0.0
	for math.Abs(_z-zk) > 1e-4 {
		zk = _z
		sinp = _z / math.Sqrt(r2+_z*_z)
		v = WGS84_R / math.Sqrt(1.0-e2*sinp*sinp)
		_z = z + v*e2*sinp
	}

	// B pos[0]=r2>1E-12?atan(z/sqrt(r2)):(r[2]>0.0?PI/2.0:-PI/2.0);
	// L pos[1]=r2>1E-12?atan2(r[1],r[0]):0.0;
	if r2 > 1e-12 {
		b = math.Atan(_z / math.Sqrt(r2))
		l = math.Atan2(y, x)
	} else {
		if z > 0 {
			b = math.Pi / 2.0
		} else {
			b = -math.Pi / 2.0
		}
		l = 0
	}

	h = math.Sqrt(r2+_z*_z) - v

	//double e2=FE_WGS84*(2.0-FE_WGS84),r2=dot(r,r,2),z,zk,v=RE_WGS84,sinp;
	//for (z=r[2],zk=0.0;fabs(z-zk)>=1E-4;) {
	//zk=z;
	//sinp=z/sqrt(r2+z*z);
	//v=RE_WGS84/sqrt(1.0-e2*sinp*sinp);
	//z=r[2]+v*e2*sinp;
	//}
	//pos[0]=r2>1E-12?atan(z/sqrt(r2)):(r[2]>0.0?PI/2.0:-PI/2.0);
	//pos[1]=r2>1E-12?atan2(r[1],r[0]):0.0;
	//pos[2]=sqrt(r2+z*z)-v;

	return

}
